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AN  APL  FUNCTION  FOR  BIVARIATE  NORMAL  PROBABILITIES 


INTRODUCTION:  Bivariate  normal  distributions  have  many 
applications  such  as  in  combat  modelling,  weapons  systems 
effectiveness  studies  and  weather  prediction  problems;  several 
other  applications  are  discussed  in  [4].  It  is  well  known  that 
probability  statements  for  a  general  bivariate  normal 
distribution  can  be  transformed  into  equivalent  statements  for  a 
standard  bivariate  normal  distribution  with  0  means  and  variances 
equal  to  1.  It  is  thus  sufficient  to  be  able  to  compute 
cumulative  probabilities  for  a  standard  bivariate  normal 
distribution  with  probability  density  function 
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Similar  to  the  case  with  the  univariate  normal  distribution,  the 
bivariate  cumulative  distribution  function  (c.d.f.) 


F(h,k) 


P[XSh, YSk] 


k  h 


f(x,y) 


dx  dy 


does  not  have  a  closed  form  solution  and  numerical  integration  is 
the  usual  approach  to  evaluating  such  integrals.  Because  of  the 
importance  of  the  distribution,  the  National  Bureau  of  Standards 


(NBS)  has  published  an  extensive  set  of  tables  for  P[X>h, Y>k] 

[5].  Another  set  of  tables  using  a  different  approximation ,  has 
been  generated  by  Oven  [6].  Other  ways  to  approximate  the 
bivariate  normal  integral  are  discussed  in  [1],  [3]  and  [7].  All 
of  these  approximations  involve  numerical  integration  and  require 
extensive  computer  programming,  rendering  them  to  be  not  readily 
suitable  for  obtaining  on  the  spot  results.  These  days,  with  the 
easy  availability  of  microcomputers,  it  would  be  useful  to  have  a 
program  to  compute  bivariate  normal  probabilities  interactively 
and  also  be  able  to  incorporate  such  a  program  within  other 
application  programs.  This  will  allow  an  analyst  to  perform 
sensitivity  studies  by  varying  the  input  parameters  in  an 
application  program  and  observing  the  effect  on  the  measures  of 
effectiveness  of  interest. 

Recently,  Wang  [8]  developed  a  new  algorithm  to  compute 
bivariate  normal  probabilities,  that  does  not  require  numerical 
integration,  relatively  easy  to  program  and  provides  quite 
accurate  results.  Investigations  by  Wang  indicate  that  the 
computed  probabilities  compare  very  well  with  those  in  the  NBS 
tables  and  the  computer  resources  needed  are  not  excessive.  The 
approach  used  by  Wang  is  to  start  with  an  approximate  contingency 
table  of  cell  probabilities  for  a  rectangular  grid  on  the  x,  y 
plane  and  then  apply  the  "iterative  proportional  fitting 
algorithm"  [3;  sec  3.5]  to  modify  the  cell  probabilities;  the 
iteration  is  continued  until  the  marginal  probabilities  in  the 
contingency  table  coincide  with  their  true  values  to  within  a 
specified  degree  of  accuracy.  Since  the  marginal  distributions 
of  the  random  variables  X  and  Y  are  univariate  normal,  the  exact 


marginal  probabilities  can  be  determined  using  computer  programs 
available  in  most  statistical  software  packages  or  by  writing  a 
subprogram  for  the  purpose.  It  should  be  noted  that,  although  a 
bivariate  normal  distribution  is  defined  over  the  entire  plane, 
for  all  practical  purposes  it  is  sufficient  to  consider  only  the 
square  subregion  [-4,4]  x  [-4,4]  since  the  probability  content 
outside  of  this  subregion  is  negligible. 

Wang's  algorithm  is  defined  by  the  following  steps: 

1.  Let  -4  ■  aQ  a1  ...  a,,  *  4  be  a  partition  of 

the  interval  [-4,4]  along  the  x  -  axis  and  -4  -  b0  bx 
...  bn  -  4,  a  partition  along  the  y  -  axis;  identical 
equal  length  partitions  for  both  axes  is  recommended  to 
reduce  computational  complexity. 

2.  Let  Ax  -  x^  -  a^_^  i  ■  1,  2,...,  m 

Ay  -  bj  -  bj.i  j  -  1,  2,...,  n 
m^  ■  a^..^  +  a^  i  ■  l,2...,m 

nj  "  bj-i  +  bj  ^  *  1,2, ... ,n 

P  -  P  (1  -  Ax^)  (1  -  Ay2  ) 

12  12 

where  is  the  correlation  coefficient  (specified) . 

3.  Let  Pi  *  p  [ai_1  <  x  <  a^l  i  ■  1,2,  —  ,m 

-p3  ■  MV1  <  Y  <  bj]  J.1,2 . n 

be  the  marginal  probability  contents  of  the  i-th  and  j-th 
subintervals  along  the  x  and  y  axes  respectively. 


4 


.  Compute  ,  the  starting  approximate  probability  of  the 

(i,  j)  th  cell  as 


p  ainj 
1>P2 


i  ■  1, 2 , . . .m;  j  *  1,2, ... ,n 


5.  The  application  o t  the  iterative  proportional  fitting 
algorithm  results  in  the  following  equations  for  the  modified 
probability  of  the  (i,j)  th  cell  after  the  k-th  iteration: 


,00 

ij 


(k-1) 


V  p(k‘1} 

X, 

j-1 

(k-1) 

P  P 

ij  j 

tC’ 

i-1 


k  -  1,3,5, . . . 


k  «  2,4,6, . . . 


6.  Continue  the  iteration  process  until  for  some  even  number  k 
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where  e  is  a  prespecified  degree  of  accuracy  with  which  the  true 
marginal  probabilities  agree  with  the  marginals  in  the 
contingency  table. 

7.  To  compute  P[X_<h,Y^k]  (or  P  [  X  >  h,  Y  >  k  J  )  sum 
the  probabilities  in  the  contingency  table  over  those  cells  for 
Which  a  <  h  (a.  >  h)  and  b,  <  k  (b.  >  k) .  In  those  cases  where  either 

i  ”  1  •  J  J 

h  +  Si  and/or  k  +  b  for  any  of  the  partition  points  a^^  and  ,  the  j 
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accuracy  of  the  approximation  can  be  improved  by  including  h 
and/or  k  as  additional  partition  points. 

An  APL  function,  called  BVN,  to  compute  bivariate  normal 
probabilities  (both  P  [  X  £  h  ,  Y  £  k  ]  and  P  [  X  >  h,  Y  >  k  ])  is 
presented  in  the  appendix.  This  function  invokes  another  APL 
function  called  NCOF  to  compute  the  marginal  univariate  normal 
probabilities.  The  BVN  function  can  be  run  on  an  IBM-PC  /  AT 
compatible  microcomputer  using  an  APL  language  system  such  as 
APL* PLUS  from  the  STSC  corporation;  the  function  can  also  be  run 
under  VSAPL  on  the  NPS  mainframe  computer.  The  function  runs 
interactively  and  calls  for  keyboard  input  of  the  desired  equal 
length  partition  size  for  the  x  and  y  axes  and  the  degree  of 
accuracy  e  in  approximating  the  marginal  cell  probabilities. 
With  only  minor  modifications,  the  function  can  be  imbedded 
within  another  APL  function  as  a  subprogram.  Computations  using 

the  BVN  function  indicate  that  with  partition  size  x  -  y  -  .2 

for  both  the  x  and  y  axes  and  e  -  .00005,  a  four  decimal  place 

accuracy  as  compared  with  the  NBS  tables  can  be  achieved.  The 

computational  time,  as  is  to  be  expected,  increases  with  a 
decrease  in  the  partition  size  x  or  y,  an  increase  in  and 
to  a  lesser  degree  a  decrease  in  e  .  With  x  »  y  ■  .2, 

c  =  .00005  and  .  1£  <.8  the  computational  time  (clock  time) 

was  between  90  and  180  seconds  on  the  Zenith  Z-248  (an  AT  type) 
microcomputer,  and  on  the  IBM  3033  mainframe  computer  these  times 
were  between  1  and  20  seconds.  The  computational  time  can  be 
reduced  considerably  (to  about  30  seconds  on  the  Z-248)  by 
choosing  x  ■  y  «  .5  but  then  only  a  three  decimal  computational 


accuracy  can  be  expected.  For  a  fixed  value  of  if  the  c.d.f. 
is  to  be  calculated  for  several  choices  of  (h,k)  ,  the  BVN 
function  needs  to  be  run  only  once;  with  a  very  minor 
modification  the  contingency  table  of  bivariate  cell 
probabilities  can  be  saved  in  a  matrix  and  all  that  is  left  is  to 
sum  the  probabilities  of  the  appropriate  cells.  If  needed,  it 
would  be  quite  straight  forward  to  generate  tables  for  various 
choices  of  and  (h,k) . 

Professor  I.  O'muircheartaigh  of  the  O.R.  department  and  I 
are  in  the  process  of  completing  a  Fortran  program  for  the 
problem  and  expect  to  submit  the  code  for  publication. 
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APPENDIX 


THIS  APPENDIX  CONTAINS  THE  LISTING  OF  TWO  APL  VARIABLES  BVNHOW  AND 
NCDFHOW  AND  TWO  APL  FUNCTIONS  BVN  AND  NCDF.  THE  TWO  HOW  VARIABLES 
PROVIDE  SHORT  DESCRIPTIONS  OF  THE  COMPUTATIONAL  SCHEMES,  THE  INPUT 
PARAMETERS  AND  THE  SYNTAX  FOR  THE  TWO  FUNCTIONS. 

BVNHOW 


THE  FUNCTION  BVN  COMPUTES  THE  C.D.F.  OF  A  8TANDARD  BIVARIATE 
NORMAL  DISTRIBUTION  WITH  CORRELATION  COEFFICIENT  p,  U8ING  AN 
ALGORITHM  PROPOSED  BY  YUCHUNG  J.  WANG  (BIOMETRIKA,  1967,  NO. 74,  185-90). 
THE  ALGORITHM  CONSISTS  OF  PARTITIONING  THE  X-Y  PLANE  INTO  RECTANGULAR 
CELL8,  AN  INITIAL  APPROXIMATION  OF  THE  CELL  PROBABILITIES  AND  AN 
ITERATIVE  SCHEME  TO  MODIFY  THE  CELL  PROBABILITIES.  THE  ITERATION 
PR0CES8  IS  TERMINATED  AS  SOON  A8  THE  MARGINAL  PROBABILITIES  COINCIDE 
WITH  THEIR  EXACT  VALUES  (THAT  ARE  UNIVARIATE  STANDARD  NORMAL 
PROBABILITIES, COMPUTABLE  USING  THE  APL  FUNCTION  NCDF)  TO  WITHIN  A 
SPECIFIED  DEGREE  OF  ACCURACY  *  .  THE  SYNTAX  FOR  THE  FUNCTION  18 

RHO  BVN  W 

WHERE  RHO  IS  THE  CORRELATION  COEFFICIENT  AND  W  -  (x,y>  IS  THE  POINT 
AT  WHICH  THE  C.D.F.  IS  TO  BE  COMPUTED.  THE  FUNCTION  WILL  CALL  FOR 
THE  INPUT  OF  THE  DESIRED  LENGTH  FOR  AN  EQUAL  PARTITIONING  OF  THE 
INTERVAL  C~4,43  ALONG  THE  X  AND  Y  AXES  AND  c  THE  DESIRED  DEGREE 
OF  ACCURACY  IN  THE  MARGINAL  PROBABILITIES.  THE  RECOMMENDED  CHOICES 
ARE  .2  FOR  THE  PARTITION  LENGTH  6x ,  AND  .00005  FOR  «  .  THE  TOTAL 
COMPUTATION  TIME,  ON  THE  ZENITH  Z-248  MICROCOMPUTER,  SHOULD  BE  BETWEEN 
90  AND  180  SECONDS  DEPENDING  ON  THE  SIZE  OF  RHO. 

NCDFHOW 


THE  FUNCTION  NCDF  COMPUTES  THE  C.D.F.  OF  A  STANDARD  NORMAL  DISTRIBUTION, 
U8ING  THE  APPROXIMATION  DEFINED  IN  EQUATION  26.2.17,  PAGE  932  IN 
THE  HANDBOOK  OF  MATHEMATICAL  FUNCTIONS,  M. ABRAMOWITZ  AND  I.A.STEGUN 
EDIT0R8,  PUBLISHED  BY  THE  NATIONAL  BUREAU  OF  STANDARDS  (REF.  C53).  THIS 
APPROXIMATION  IS  ACCURATE  TO  ATLEAST  TO  7  DECIMAL  PLACES.  THE  SYNTAX 
FOR  THE  FUNCTION  ISl  NCDF  Z  WHERE  Z  IS  AN  INCREASING  ARRAY 
OF  NUMBERS  FOR  WHICH  THE  C.D.F.  IS  TO  BE  COMPUTED.  THIS  FUNCTION  IS 
INVOKED  BY  THE  BVN  FUNCTION  TO  COMPUTE  MARGINAL  PROBABILITTES. 

-9- 


9  RHO  BVN  W|X|Y|RBAR|R|A|M|B«E|Q|Df  IjJ|K|L1|L2|A1)A2|M1|M2|P1|P2 
C13  '  * 

C23  '  * 

C33  ft THIS  FUNCTION  COMPUTES  THE  CUMULATIVE  PROBABILITIES  OF  A  STANDARD 

C43  ftBIVARIATE  NORMAL  DISTRIBUTION  (  m«ant  0  and  st.davs  1)  WITH 

CS3  ft CORRELATION  COEFFICIENT  RHO.  THE  SYNTAX  FOR  THE  FUNCTION  ISi 

C63  ft RHO  BVN  W  ,  WHERE  W-<xfy).  THE  FUNCTION  FIRST  COMPUTES  A 

C73  ft CONTINGENCY  TABLE  OF  CELL  PROBABILITIES  OVER  A  USER  DEFINED  PARTITION 

C83  ft OF  THE  X  AND  Y  AXES  AND  THEN  CUMULATES  THE  PR0BAILITIE8  IN  THE 

C9]  ft APPROPRIATE  CELLS.  THE  USER  IS  PROMPTED  TO  INPUT  THE  LENGTH  OF  THE 

C 103  ftGUBINTERVALS  IN  THE  PARTITION  OF  (“4, 4)  («.g. ,  .05)  AND  THE  DESIRED 

C113  ft  ACCURACY  <«.g.,  .0005  FOR  A  3-DECIMAL  ACCURACY).  THE  FUNCTION 

C 123  ft OUTPUTS  BOTH  PrCX<x,Y<y3  AND  Pr  C  X  >  x  fY>y3. 

C 133  '  INPUT  THE  DESIRED  PARTITION  SUBINTERVAL  LENGTH  FOR  X  AND  Y  AXES' 

C 143  K*0 

C1S3  '  ' 

C 163  'INPUT  THE  DESIRED  ACCURACY  OF  COMPUTATIONS' 

C173  E+0 

C 183  X+WC 1 3 
C 193  Y+WC23 

C203  ♦END1x\(XS“4)vY1“4 

C213  ♦END2xt <X*4)~Y*4 
C223  +END3XV (Xl4)vY^4 
C233  X*L/4f  X 
C243  Y«-L/4,Y 

C253  RBAR«-RHOxl— <K*2) *12 
C263  R*RBAR* ( 1 — RBAR*2) 

C273  A^"4+OfKXT.0*K 

C283  AIM  <X>A)/A)  ,X,  <X<A)/A 

C293  A2M  <Y>A) /A) fY, <Y<A) /A 

C303  Ml«-<  Cl*Al)+“l*Al)+2 

C313  M2M  <14A2)+”l+A2)+2 

C323  Ll«-pMl 

C333  L2+PM2 

C343  B^*RxM1o.xM2 

C353  PI*  (NCDF  UAD-NCDF  “14A1 

C363  P2MNCDF  UA2J-NCDF  “1+A2 

C373  REPEAT  1  Q*P1**/B 

C383  IK< (Ll,Ll)pQ)x( (LI ,Ll)pl ,LlpO) 

C393  B+D+.xB 
C403  Q+P2+W  B 


[413  D*<  (L2,L2)/>Q)x( (L2f L2)pi f L2p0) 

[423  B+B+.xD 

[433  +REPEATx\ <  <+/ < I P1-+/B) >E>  >0>  v  <+/ < I P2-+/B) >E) >0 
[443  I«h-/X>A1 

[453  J«*+/Y>A2 

*  [463  'Pr  t  X  <  'f(fX>,'  ,  Y  <  *  f  CfY) ,  *  3  -  ' , f+/+/ <1 f J)+B 

[473  'Pr  [  X  >  ' ,  <»X> , '  ,  Y  >  ' ,  <»Y>  f '  3  -  ',»■»■/+/  (I  fJ)*B 

*  [483  +0 

C493  ENDli '  Pr  [  X  <  ',C»X)f'  f  Y  <  ',(»Y>,'  3  -  0' 

[303  '  Pr  [  X  >  ',<¥X>f'  ,  Y  >  '  f  <¥Y>  , '  3  -  1* 

[313  *K> 

[323  END2I '  Pr  [  X  <  ',<¥X)f'  ,  Y  <  ',(»Y),‘  3  -  1 ' 

[533  '  Pr  [  X  >  ',(¥X)f'  ,  Y  >  ',<¥Y>f'  3-0' 

[543  +0 

[533  END3I '  Pr  [  X  <  'f<¥X>,'  ,  Y  <  ',(*Y)f'  3  »  ' ,* (NCDF  X)xNCDP  Y 

[563  '  Pr  [  X  >  ',CfX)t'  f  Y  >  ' f <¥Y> , '  3  -  'f¥Cl-NCDF  X)x(l-NCDF  Y> 

[373 

v 
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7  RH4C0F  Z|B|p|Z|T|N|M|P 

£13  A THIS  FUNCTION  COMPUTES  THE  C.D.F.  Pr£ZSz30FA  STANDARD 
£23  A NORMAL  DISTRIBUTION  USING  THE  APPROXIMATION  IN  EQUATION  26.2.17, 
C33  ft PAGE  932  IN  THE  HANDBOOK  OF  MATHEMATICAL  FUNCTIONS,  EDITED  BY 
£43  AABRAMONITZ  AND  STAT6UN,  NATIONAL  BUREAU  OF  STANDARDS. 

CS3  Z+,Z 
£63  Z+ZE4Z3 
£73  N*p<Z<0)/Z 

£83  M*p<Z20)/Z 

£93  Z«-IZ 
£103  p+0. 2316419 

£113  T«"H+p*Z 

£123  *♦<*—<  Z*2) +2) ♦ <o2>  *0. 5 

£133  B*  0.31930133  “0.336363782  1.781477937  *1.821235978  1.330274429 
£143  P«*l— zx  (To  .  *x3>  +.  xB 

E 133  R*  < 1— NfP) , <-M>  tP 
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